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From  the  principle  of  minimum  complementary  energy  which  is 
extended  to  take  into  account  the  possible  discontinuity  in  stresses  or 
displacements  along  the  interelement  boundaries,  the  finite  element 
methods  by  both  the  assumed  stress  hybrid  model  and  the  equilibrium 
model  can  be  formulated.  This  formulation  enables  the  development 
of  a  rational  method  for  consistent  lumping  of  body  and  surface  forces 
and  the  establishment  of  a  criterion  for  the  kinematical  instability  of 
a  system.  For  the  hybrid  model,  there  is  an  optimum  choice  of  the 
number  of  stress  modes  for  a  given  boundary  displacement  approx¬ 
imation.  Example  calculations  for  plate-bending  problems  are  included 
to  substantiate  the  analytical  predictions.  Examples  are  also  carried 
out  to  illustrate  the  convenience  in  taking  transverse  shear  effects 
into  account  by  the  hybrid  method. 
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SECTION  I 
INTRODUCTION 


When  the  assumed  stress  hybrid  model  was  suggested  by  Pian  (References  1  and  2)  to 
derive  the  element  stiffness  matrices,  the  complementary  energy  of  only  the  individual 
elements  was  considered.  Later  on,  a  variational  principle  was  developed  (Reference  3)  to 
treat  the  entire  assembled  elements  as  a  whole.  It  was  then  possible  to  derive  the  loading 
matrix  in  a  consistent  manner  and  to  reveal  some  restrictions  in  the  formulation  of  the 
assumed  stress  hybrid  method. 

The  present  paper  is  to  derive  a  variational  principle  which  can,  in  fact,  cover  both  the 
assumed  stress  hybrid  model  and  the  equilibrium  model  by  Fraeijs  de  Veubeke  (Reference  4 
and  5);  hence  some  of  the  characteristics  are  common  to  both  methods.  The  formulation  of 
the  hybrid  model  is  then  extended  to  plate-bending  problems  including  both  thermal  effect 
and  transverse  shear  effect.  Example  calculations  are  made  to  evaluate  this  method. 
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SECTION  n 

PRINCIPLE  OF  MINIMUM  COMPLEMENTARY  ENERGY- 
EQUILIBRIUM  MODEL  AND  ASSUMED  STRESS  HYBRID  MODEL 


Both  the  equilibrium  model  and  the  assumed  stress  hybrid  model  for  finite  element 
analysis  can  be  derived  from  the  principle  of  minimum  complementary  energy,  for  which  the 
functional  to  be  varied  is 


In  this  expression 


a , 


(I) 


^ij  =  stress  tensor 

Cijki=  elastic  compliance  tensor 

V  =  volume 

S  =  surface 

T.  =  surface  traction 

u.  =  prescribed  boundary  displacement 

®u  ~  portion  of  S  over  which  the  boundary 
displacements  u.  are  prescribed 


The  stress  tensor  <x  satisfies  the  equilibrium  conditions 

cr. .  .  +  F,  =  0  (2 

•  i  ,  J  i  ' 

and  is  compatible  with  the  prescribed  boundary  tractions  over  the  boundary  S  .  F.  is  the 
prescribed  body  force.  ! 


In  applying  the  finite  element  method,  the  assumed  stress  field  need  not  be  continuous 
across  the  interelement  boundaries,  but  equilibrium  must  be  maintained  for  the  surface 
tractions  T.,  defined  by  T.  =  cr..  v.,  where  v  are  the  components  of  the  unit  vector  normal 
to  the  boundary.  Let  two  neighboring  elements  “a”  and  “b”  be  isolated  (Figure  1),  and  con¬ 
sider  the  boundary  traction  components  T.(a)  (s)  and  T.(b>  (s)  (i  =  1,2,3)  over  the  respective 
sides  of  the  common  boundary  AB.  (For  simplicity,  a  plane  stress  example  is  used  for 
illustration.) 
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Figure  1.  Equilibrium  of  Boundary  Tractions  Along  Interelement  Boundary 


The  equilibrium  conditions  at  the  interelement  boundary  are  given  by 

TU)(s)  +  T(.b)(  s )  5  0  (  i  =1,2,3)  (3) 

i  i 

Equation  3  may  be  considered  as  conditions  of  constraint  and  can  be  introduced  by  including 
Lagrange  multiplier  terms 


or 


f  X  (s  )  |Va)(  s)  +  T(b)  {  s)  ]  ds 
JAB  1  L  1 


/ 


T: 


X:  T : 


(  4) 


in  the  complementary  energy  functional  to  be  varied.  The  Lagrange  multipliers  X^,  which  are 
functions  of  the  surface  coordinates,  are  to  be  treated  as  additional  variables. 


When  all  the  interelement  boundaries  have  been  considered,  the  complementary  energy 
functional  may  be  written  as: 


(5) 
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where  S  is  the  boundary  of  V  where  the  displacements  are  prescribed  and  s  is  the 
n  n  n 

interelement  boundary  of  V  .  By  taking  the  variation  of  7 r  with  respect  to  a. .  and  X . ,  it  is 

n  i]  i* 

easily  shown  that  the  Lagrange  multipliers  \.  are  equal  to  u.,  the  displacements  along  the 
interelement  boundary 


?  (/ 

vn 

-/ 


~  c-)ki  %  \i dv 


ijl 


Ti  “i  dS  7 

Su. 


T.  u.  dS 


(6) 


In  the  finite  element  solution,  the  assumed  approximate  functions  for  the  stresses  <x.. 
are  divided  into  two  parts.  The  first  part,  which  consists  of  a  finite  number  of  parameters  j3  , 
should  satisfy  the  homogeneous  equations  of  equilibrium,  while  the  second  part  is  a  particular 
solution  of  the  equations  of  equilibrium  with  the  prescribed  body  forces.  In  matrix  form,  the 
stresses  <r  are  expressed  as: 

<r  =  P  £  +  Pp  $F  (7) 

where  $  is  unknown  and  Pp  is  determined  from  the  particular  solution.  For  elements 
which  contain  boundaries  with  prescribed  surface  tractions,  some  of  the  /3’s  in  the  first  term 
will  also  be  prescribed.  In  this  case,  all  the  prescribed  /3  ’s  are  put  in  the  second  term. 


Both  the  equilibrium  and  the  hybrid  methods  can  be  derived  from  Equation  6,  but  the 
distinction  lies  in  the  treatment  of  the  surface  integrals  along  the  interelement  boundaries. 
For  the  equilibrium  model,  the  surface  tractions  T.  along  each  boundary  are  represented 
uniquely  by  the  generalized  loads  Q  pertaining  to  the  boundary  and  can  be  written  as 

T  =  <f>  Q  (8) 

For  example,  when  Q  are  the  values  of  T.  at  a  finite  number  of  boundary  points,  <f>  represents 
the  corresponding  interpolation  functions.  Since  the  surface  tractions  are  also  related  to  the 
assumed  stress  distribution,  they  can  also  be  expressed  as 

T  -  R/9  -b  rfj3f 

Thus,  there  is  a  unique  relation  between  Q  and  /3  of  the  form: 

o  .  gt  $  +  gtf  0f 

The  corresponding  element  generalized  displacements  <1  are  defined  by: 


(  9) 


(10) 


T 

Q  q 


■/ 


T.  u.  dS 


>7  *r 


«B  dS 


(ID 
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or  r  i  T 

q  -  J  4>  Ug  dS  (12) 

Sn 

Thus,  the  generalized  displacement  is  a  weighted  integral  of  the  boundary  displacements,  and 
to  maintain  the  same  q  for  two  neighboring  elements  does  not  guarantee  compatibility  along 
the  entire  boundary.  It  is  seen  that,  in  an  equilibrium  model,  each  generalized  displacement 
is  common  to  only  two  neighboring  elements. 

For  the  hybrid  method,  the  approximate  displacements  along  the  interelement  boundaries 
are  represented  by  interpolation  functions  and  the  generalized  displacements  q  at  a  finite 
number  of  boundary  nodes: 

uB  =  L  q  (13) 

Since  the  interpolating  functions  L  are  applied  to  the  individual  boundary  segments,  they  are 
relatively  easy  to  construct  so  that  interelement  compatibility  is  maintained.  Unlike  the 
equilibrium  model,  the  generalized  displacements  q  may  be  referred  to  corner  nodes  where 
more  than  two  elements  are  connected. 


The  corresponding  generalized  nodal  forces  are  defined  again  by: 

qT  Q  =  /  T.  Uj  dS  (14) 

"g 

or  n 

Q  -  f  L  T  dS  (15) 

S  n 

Thus,  the  generalized  nodal  force  is  a  weighted  integral  of  the  boundary  tractions,  and  to 
maintain  equilibrium  at  a  node  does  not  guarantee  equilibrium  over  the  entire  boundary. 
Another  significant  difference  from  the  equilibrium  model  is  that  the  stresses  and  the  boundary 
displacements  are  independently  assumed.  Hence,  the  number  of  parameters  in  /3  and  the 
number  of  generalized  displacements  q  can  be  chosen  independently  in  the  hybrid  model. 


Finally,  since  the  hybrid  method  relies  on  assumed  boundary  displacements,  the  pre¬ 
scribed  boundary  stresses  no  longer  constitute  a  restrained  boundary  condition.  In  this  case 
the  functional  tt  can  be  written  more  conveniently  as 


7T  ; 


—  c 
a  ijk  i? 


cr  . 
'  J 


V 


dV 


T;  U:  dS  + 


I 


T.  u  .  dS 


(16) 
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where 

a  v,  =  s„  +  Strn  +  su„ 

is  the  entire  boundary  of  V  ,  and  u.  ="u.  on  S  . 

n  i  i  u 

n 

Substitution  of  Equations  7,  10,  and  11  into  Equation  6  and  substitution  of  Equations  7,  9, 
and  13  into  Equation  16  both  yield  equations  of  the  form 

"•=  I(^-/9TH0+  /3THF£F-/3TSq 

+  s  q  +  B„)  (IT) 

where 

H  -  f  PTC  P  dV 
J  vn 

H F  =  L  pTC  PF  dV  (18) 

vn 

Bo  "±£j  Pf  CPF  dV  @F 

and  C  is  the  elastic  compliance  matrix.  For  the  equilibrium  model,  G  and  are  defined  by 
Equation  10.  For  the  hybrid  model  they  are  given  by: 

S  =  /  R*  L  dS  ;  G  =  f  Rc  L  dS  (19) 

X  F  \ 

T 

The  vector  S  is  also  different  for  the  two  models 


S  =  )9p  Gp  (equilibrium  model) 

S  =  -  0F  ®F  +/  T  t_  dS  (hybrid  model) 
S°n 

where  T  are  the  prescribed  boundary  tractions. 


(20) 


The  stationary  conditions  of  the  functional  given  by  Equation  17  with  respect  to  variations 
of  /3  and  q  then  yield 

H  @  +-  Hp  £p  -  G  q  =  0  (21) 


and 

I  (  £T  G  -  ST  )8  q  =  0 

n 


(22) 
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By  solving  for  f3  from  Equation  21  and  substituting  back  into  Equation  17,  we  can  express  the 
functional  tt  in  terms  of  the  generalized  displacements  q  only,  i.e., 

TT  :  -  5;  (1  qTt  q  -0Tq  +  cn)  (23) 

n  2 

where  k  and  Q  are,  respectively,  the  element  stiffness  matrix  and  equivalent  nodal  forces 
defined  by 

k  =  GT  H_l  G  (24) 

Q  =  GT  H"'  Hp  )9F+  S  (25) 

and 

T  t  _ 

Cn  -  “2”^p  Hp  H  1  Hp  /?p  -  Bn  -  Constant  (26) 

Knowing  k  and  Q  for  all  elements,  it  is  then  a  routine  matter  to  set  up  the  assembled 
matrix  equation  for  the  problem. 

Inspection  of  Equation  22  reveals  that  if  N,  the  total  number  of  assumed  stress  modes  ft 
of  all  elements  is  smaller  than  M,  the  total  number  of  unknown  displacements  q  ,  there  will 
be,  in  general,  no  solution  for  the  /3’s.  It  can  be  seen  that  the  above  situation  corresponds  to 
the  appearance  of  kinematical  deformation  modes  discussed  by  Fraeijs  de  Veubeke  (Ref¬ 
erence  5),  who  indeed  had  pointed  out  that  for  a  single  element,  if  n  is  the  total  number  of 
stress  modes  and  m,  the  number  of  generalized  displacements  of  which  r  must  be  constrained, 
r  being  the  rigid  body  degrees  of  freedom,  then  for  such  element  to  be  kinematically  stable, 
n  must  be  larger  than  (m  -  r).  It  is  obvious  that  such  requirements  on  the  minimum  number 
of  stress  modes  can  always  be  met  by  the  hybrid  model  because  the  stress  modes  and 
boundary  displacements  are  assumed  independently,  but  it  is  not  so  for  the  equilibrium  model. 
Furthermore,  since  for  the  hybrid  model,  there  will  always  be  more  than  two  elements 
meeting  at  an  internal  node,  the  assembled  structure  will  be  kinematically  stable  if  the  above 
criterion  for  each  individual  element  is  satisfied.  For  the  equilibrium  model,  however,  an 
assembled  structure  may  be  unstable  even  if  all  the  individual  elements  satisfy  the  criterion 
for  kinematical  stability  (Reference  5). 
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SECTION  HI 

COMPLEMENTARY  ENERGY  PRINCIPLE  FOR  PLATE  THEORY 
INCLUDING  THERMAL  AND  TRANSVERSE  SHEAR  EFFECTS 


The  thermal  effect  is  usually  included  in  the  variational  principle  as  an  initial  strain 
problem  (Reference  6).  Let  be  the  initial  strain  and  e  ..  the  strain  due  to  total  de¬ 

formation,  then  the  stress  is  given  by 


°"i  j  E  i  j  ^  €k l  €oki  * 


(27) 


It  is  convenient  to  define 


ij  =  Eijk^  €ki  =  stress  associated  with  the  total  deformation 


°o  i  j  E  i  j  k  i  €o  k  i 


=  stress  associated  with  the  initial  strain 


Then 


cr  =  cr  -  cr  . . 
'j  «j  01 J 


(28) 


(29) 


For  the  plate-bending  problem,  it  is  more  convenient  to  use  cr/.  as  the  primary  variable 
because  the  in-plane  stress  cr^^g  varies  linearly  across  the  thickness;  hence  can  be  defined 
uniquely  by  the  stress  resultants  and  stress  couples  M ,  i.e., 

,  Kb 

aa/3  1 


1 2  M 


a(3 


h  h3 


(30) 


Here  the  Greek  index  takes  the  value  of  1  and  2  and  the  coordinate  normal  to  the  mid-plane 
is  Xg.  (Figure  2) 


For  the  plate  theory  including  transverse  shear  effect,  the  development  in  the  following 
is  based  on  the  treatment  of  Reissner  (Reference  7). 


Substituting  Equations  29  and  30  into  the  equilibrium  equations  in  three-dimensional 
elasticity 

%,i  =0  13,1 

for  i  =  1,2  and  integrating  through  the  thickness  direction  we  obtain  the  expressions  for  the 
transverse  shear  stress  which  satisfy  the  shear  free  conditions  on  both  faces  of  the  plate  is 


03 


(32) 
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Figure  2.  Sign  Convention  for  Plate-Bending  Problem 


where 


Ja 


h/4 

\  h  2  /  V  u,  c 


-n/2 


oa/3,  /3  d*3 


-X 


°oa/3,/3 


dx. 


Qa  =  M' 


=  transverse  shear 


(33) 


(3  4) 


'a/3,/3 

The  substitution  of  Equation  32  into  Equation  31  for  i  =  3,  and  integration  along  Xy  yields 

with  a„„(-h/2)  =  0.  By  introducing  the  condition  a-  (  +  h/2)  =  p,  we  obtain 

OO 


(35) 


h/2 

Qa,ct=  _p  +/h/2Ba,a 


dx. 


(36) 


It  is  obvious  that  the  terms  involve  Ba  in  Equation  32  and  Bg  a  in  Equations  35  and  36 
are  the  stresses  due  to  thermal  effect  in  the  present  plate  theory 


^oas  =  Ba 


(37) 
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The  complementary  energy  functional  can  be  expressed  in  terms  of  either  cr. .  —  ' 


i.e., 


..  or  cr... 
ij  ij 


h4 

'«■/  I. 


'A  2  Cijkj?  °ij  ‘V'Sij  °ij  >1*  <l«. 


-•f/S 

-h/2 


cr. .  z/.  u.  ds  dx, 
ij  i  i  3 


(38) 


or 


-h/2 


%i  dAd‘=  -/ 


u 


ct-.  !  i/.  u.  ds  dx_ 
su  ij  j  i  3 


(39) 


By  limiting  our  discussion  to  bending  problems  only  and  by  expressing  the  stress  due  to 
deformation  in  the  same  form  as  that  of  Reference  7 


t 

/T*  /s  - 

12  Mq/3  x. 

°a/3 

h5 

/ 

3Qa  r , 

a  - 
az 

2  h  [" 

/ 

3  r  2 

a33 

T  p  IT 

llL)2' 

h  '  ■ 
2x„ 


(40) 


+ 


-  ] 


Substituting  into  Equation  39  and  carrying  out  the  integration  with  respect  to  x„  we  obtain 

O 


the  complementary  energy  functional  in  terms  of  M^and  Q^, 


I  I  O  2 

h  C 


'A  h3  L''ap\e'"aP”\er5  "  ua3/33  Qa°/3 


+  4. 


5Ca/Q33Ma/9  P  +  -i“Ca333  Qa  p]dA 

~I  [Ma/3  v/3^a  ~  Qa  va  "  ]  ds 


(41) 


where  </>Q  is  the  rotations  such  that 
and  w  is  the  normal  displacement  of  the  plate. 


We  can  also  express  the  complementary  energy  functional  in  terms  of  MQp  and  Qa<  By 
defining  the  stress  couples  and  transverse  shears  as 

_h/2  Wz 


M 


a/3 


■/, 


IV  *3  llx3 


;  “a  *  /„.  °a 


-"A 


3  dX3 


(42) 
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and  by  writing 


we  can  easily  identify  that 


a/3  =  Ma/3 


Qa  Qoa 


r 

'oa  /3  ''  x3  aoa(3  dx3 

-  /9 


h/z 

-  /2 


dx 

o  as  * 


and  that  Qa  and  Ma£  satisfy  the  conventional  plate  equilibrium  equations.  In  terms  of  Ma/g 
and  Qa,  the  functional  7 is 

Tc  [Ca/3^V^5hZcas^Qa  Qj3 

+  -T  ca/3  3s  Ma£  p  +  ^as/^  °a  p]  dA 

-  Ma/3  M°\e+f  c«/3,0o  VJd4 

■/  lMa,8  ‘  Va  “  1  ds  (45' 

Sa 

In  the  following  special  cases,  the  expressions  for  the  complementary  energy  can  be 
simplified.  For  example: 

1.  Isotropic  material  with  thermal  strain 

77"c  s  ~  J  [(Mx  +My)%2(l^)(M2xy-MxMy) 

A  Eh 

+•  —  ff  (Q  x  +  Qy  )  -  vp  (  Mx  +M  y  )jd A 

-f  (Mx+My)  dA 

A  Eh 


1M»  *0  +  Mns  ^  *’ 
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where 


h/2 


«o=  / 

-h/2 


zE  ct  AT  (I  -Zv) 


dz 


I  —  i/2 

2.  Conventional  plate  theory  with  thermal  strain  but  with  transverse  shear  effect 
-  if  K  +Myf  +2  (l+!/)(Mx2y  -  MxMy  )j  dA 


~  f-^l^~v)MQ  (  M  x  +  My  )  d  A 
~/  (  Mn  w,  n  -f- Mns  w,  s  -  Qn  w  ) 


ds 


3.  Sandwich  plate  (with  facing  material  thickness  t^«h)  with  thermal  strain  effect 


7 r  = 
c 


=  “  /  2  g  [(Mx4MyH2(H-y)(Mxy-M]<My  H 

A  ' 

r  2 

-  - sr  (l-i/)  M0  (Mx+Mv)  dA 

•'a  Etf  h2  o  x  y 

-/  (  M  n  +  Mns  -  Qn  w)  ds 


n  Etf  1/n2 


(q:+q:)  dA 


2  Gh  x  y 


>] 


(47) 


(48) 


(49) 


where 


E  ht, 


M 


n  =  - —  a  (AT,.  -AT a  ) 

0  2(1+1/)  u  * 


The  subscripts  u  and  Jl  indicate  the  upper  and  lower  facing,  respectively. 

It  can  be  seen  that  when  Equations  45,  46,  48  or  49  are  used  in  the  formulation  of  the 
assumed  stress  finite  element  methods,  an  equation  in  the  form  of  Equation  17  will  result. 

rr 

Here  the  term  P  HF£p  will  arise  not  only  from  the  particular  solution  for  stresses  but 
also  from  the  terms  involving  Ma^gp,  Qap,  MapM.Q-^Q  and  QQQ0fi-  Noting  that  the  matrix 
H  „  appears  in  the  generalized  nodal  forces  (see  Equation  25)  we  can  conclude  that  the  initial 

r 

strains  will  lead  to  equivalent  nodal  forces  and  the  distributed  load  p  will  make  further  con¬ 
tribution  to  the  generalized  forces  when  the  “moderately  thick”  plate  theory  is  employed. 
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SECTION  IV 

BOUNDS  TO  DIRECT  INFLUENCE  COEFFICIENTS 


Finite-element  analyses  by  compatible  and  equilibrium  models  will  yield  bounds  on  the 
total  strain  energy  content.  Since  the  direct  flexibility  influence  coefficient  (i.e.,  the  gener¬ 
alized  displacement  due  to  the  corresponding  generalized  force  of  unit  magnitude)  is  equal 
to  the  strain  energy,  its  upper  and  lower  bounds  can  also  be  established  by  the  equilibrium 
and  compatible  models,  respectively.  This  idea  has  been  explored  by  Fraeijs  de  Veubeke 
(Reference  4)*  in  dual  finite  element  analyses  of  different  types  of  structures. 

The  direct  influence  coefficient  obtained  by  the  finite  element  hybrid  method  may  be 
either  an  upper  or  a  lower  bound,  since  the  method  is  not  based  on  a  minimum  principle. 
However,  it  can  be  shown  that  the  direct  influence  coefficient  by  the  hybrid  model  is  always 
bounded  by  that  of  a  compatible  model  using  the  same  type  of  interelement  boundary  dis¬ 
placements  and  that  of  an  equilibrium  model  using  the  same  type  of  interior  stresses. 

We  observe  that  within  each  element,  the  compatible  model  satisfies  the  compatibility 
conditions,  while  the  hybrid  model  satisfies  internal  equilibrium.  Therefore,  a  compatible 
model  is  more  rigid  than  the  hybrid  model  if  the  same  set  of  generalized  coordinates  q  and 
the  same  interpolating  functions  for  the  interelement  boundary  displacements  are  chosen. 
In  that  case,  the  direct  influence  coefficient  obtained  by  the  hybrid  model  will  be  larger  than 
that  obtained  by  the  compatible  model.  For  the  compatible  model,  we  can  introduce  additional 
displacement  modes  but  keep  the  same  set  of  generalized  coordinates  by  the  static  conden¬ 
sation  process  (Reference  8).  The  introduction  of  such  additional  modes  is  to  improve  the 
satisfaction  of  the  internal  equilibrium.  It  is  understood  that  the  internal  compatibility  con¬ 
dition  is  always  satisfied.  On  the  other  hand,  the  use  of  more  and  more  stress  modes  /3  in  the 
hybrid  model  improves  the  satisfaction  of  internal  compatibility,  while  the  internal  equilibrium 
conditions  are  already  satisfied.  It  is  obvious  that  the  two  methods  will  yield  exactly  the  same 
result  if  the  assumed  modes  in  both  cases  approach  infinity.  We  can  also  interpret  the  use 
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of  more  displacement  modes  in  the  compatible  model  as  a  relaxation  of  artificial  constraints. 
However,  the  use  of  more  stress  modes  in  the  hybrid  model  would  produce  a  structure  more 
rigid  than  one  with  fewer  stress  modes. 

For  the  equilibrium  model,  when  a  set  of  assumed  stress  modes  is  chosen,  there  is  only 
one  choice  of  generalized  displacements.  If  the  same  stress  modes  are  used  in  a  hybrid  model, 
the  structure  will  be  more  rigid  than  the  one  analyzed  by  the  corresponding  equilibrium  model 
because  of  the  imposed  compatibility  conditions  along  the  interelement  boundaries.  The  direct 
influence  coefficient  of  the  hybrid  method  will  be  lower  than  that  of  the  equilibrium  method 
using  the  same  stress  modes.  Table  1  gives  a  summary  of  the  conditions  fulfilled  by  the 
various  finite  element  schemes.  They  are  arranged  according  to  increasing  rigidity  of  the 
different  models. 

Since  the  finite  element  compatible  model  always  yields  a  lower  bound  of  the  direct 
influence  coefficient,  it  appears  that  the  direct  influence  coefficient  obtained  by  using  more 
internal  nodes  is  always  superior  to  that  obtained  by  using  fewer  or  no  internal  nodes.  No 
such  conclusion  can  be  made  concerning  the  number  of  stress  modes  in  the  hybrid  model.  In 
fact,  if  the  direct  influence  coefficient  obtained  by  using  i  stress  modes  in  that  method  is 
already  a  lower  bound,  then  the  use  of  j  stress  modes  with  j  >  i  can  only  make  the  solution 
more  inaccurate.  Thus,  there  is  an  optimum  number  of  assumed  stress  modes  in  the  hybrid 
method. 

There  does  not  seem  to  exist  any  method  that  one  can  use  to  predetermine  the  optimum 
number  of  assumed  modes.  The  situation  will  be  different  for  different  problems  and  for 
solutions  at  different  locations  of  the  structure.  However,  one  can  see  immediately  that  the 
resulting  error  of  the  approximate  solution  may  arise  from  both  the  error  of  the  assumed 
stress  field  and  that  of  the  assumed  boundary  displacement.  One  can  thus  conclude  that  the 
two  independent  assumed  functions  should  be  comparable,  i.  e.,  that  they  should  provide  errors 
of  the  same  order  of  magnitude.  For  example,  in  the  plate-bending  problem,  if  the  assumed 
moment  distribution  in  the  interior  of  each  element  is  a  quadratic  function  of  the  space 
coordinates,  the  transverse  shear  distribution  will  be  linear,  then  the  order  of  magnitude  of 
error  of  the  boundary  integral  /Q  wds  will  be  comparable  with  that  of  /Mnw  ^ds  when  w  n  is 
assumed  to  vary  linearly  along  the  boundary  (Reference  3).  The  order  of  magnitude  of  the 
error  cannot  be  reduced  by  including  many  assumed  stress  modes  for  the  interior  of  an 
element  without  improving  the  boundary  displacement  function.  Indeed,  in  the  case  of  the 
hybrid  model,  the  use  of  additional  stress  modes  may  actually  decrease  the  accuracy  of  the 
solution.  This  argument  of  comparable  errors  in  the  assumed  interpolating  functions  is 
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also  applicable  to  the  compatible  model.  The  use  of  many  additional  internal  modes  for  a 
given  boundary  displacement  function  again  cannot  reduce  the  order  of  magnitude  of  the 
error. 


SECTION  V 

RESULTS  OF  FINITE -ELEMENT  PLATE-BENDING  ANALYSES 

To  illustrate  the  various  points  on  the  assumed  stress  methods  brought  up  in  this  paper, 
several  numerical  results  by  finite-element  methods  are  presented  in  this  section, 

1.  The  central  deflection  of  a  square  plate  with  clamped  edges  under  central  loading 
has  been  analyzed  by  different  finite -element  methods  using  right-triangular  elements.  For 
comparison,  the  compatible  model  by  Clough  and  Tocher  (Reference  9),  the  equilibrium  model 
by  Fraeijs  de  Veubeke  and  Sander  (Reference  10)  and  the  hybrid  model  by  Severn  and  Taylor 
(Reference  11)  are  discussed.  The  results  are  plotted  in  Figure  3  in  terms  of  the  percentage 
error  versus  the  number  of  meshes  in  half  of  the  plate  width.  The  results  clearly  indicate 
the  upper  and  the  lower  bound  solutions,  respectively,  by  the  equilibrium  and  compatible 
models.  It  is  seen  that  since  in  both  the  compatible  model  and  the  hybrid  model  the  inter¬ 
element  displacement  functions  are  cubic  in  w  and  linear  in  w  ,  the  latter  model,  according 
to  our  previous  conjecture,  yields  a  more  flexible  structure.  On  the  other  hand,  since  in  the 
hybrid  model  the  stress  approximation  involves  quadratic  terms  in  the  moment  while  the 
equilibrium  model  relies  on  a  less  accurate  linear  distribution,  the  hybrid  model  leads  to  a 
comparatively  more  rigid  structure.  It  should  also  be  observed  that  the  hybrid  method  pro¬ 
vides  a  lower  bound  in  this  case  and  hence  provides  more  accurate  solutions  than  the  corre¬ 
sponding  compatible  model. 

2.  The  second  example  shown  in  Figures  4  and  5  is  the  central  deflection  of  a  simply- 
supported  square  plate  under  central  load.  Figure  4  shows  the  results  obtained  by  using 
12  degrees  of  freedom  rectangular  elements.  However,  in  this  case,  an  equilibrium  model 
will  result  in  an  element  with  kinematical  deformation  modes.  Hence,  only  the  results  of  the 
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Figure  5.  Center  Deflection  of  Simply-Supported  Square  Plate  Under  Center  Load 
{Assumed  Stress  Hybrid  Method  -  Rectangular  Elements) 
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compatible  and  hybrid  models  are  presented  for  comparison.  The  formulation  of  the  hybrid 
method  for  a  rectangular  plate  element  employed  different  moment  approximations  in  the 
interior  of  each  element.  The  assumed  boundary  displacements  are  cubic  in  w  and  linear 
in  w  .  The  compatible  model  is  the  one  used  by  Deak  and  Pian  (Reference  12),  hence  involves 
the  same  boundary  displacement  approximation  as  the  hybrid  model. 

It  is  clearly  indicated  here  that  the  hybrid  model  may  yield  either  the  upper  or  the  lower 
bound  for  the  direct  flexibility  influence  coefficient,  and  with  the  use  of  a  sufficiently  large 
number  of  stress  modes,  it  will  yield  a  lower  bound  solution.  In  the  present  case,  the  use  of  a 
quadratic  moment  distribution  yields  the  most  accurate  results. 

The  hybrid  solutions  include  two  different  approaches:  in  one  case  the  assumed  stress 
distribution  satisfies  the  condition  of  vanishing  normal  moments  along  the  simply- 
supported  edge;  in  the  other,  the  prescribed  boundary  stress  condition  is  not  considered  as 
a  restrained  condition.  As  indicated  in  the  results,  the  first  of  these  two  approaches  provides 
more  accurate  solutions,  but  the  improvement  becomes  diminishingly  small  when  a  large 
number  of  elements  is  used  in  the  analysis. 

Figure  5  compares  the  results  for  the  same  problem  obtained  by  using  different  com¬ 
binations  of  interior  stress  approximations  and  boundary  displacement  approximations.  In 
addition  to  the  12  degrees  of  freedom  elements, the  16  degrees  of  freedom  elements  are  used. 

For  the  latter,  both  w  and  w  vary  as  cubic  functions  along  the  edges.  Also  in  these  sets  of 

»n 

solutions,  the  assumed  stress  modes  do  not  satisfy  the  prescribed  boundary  condition  of 
vanishing  normal  moment.  This  figure  shows  clearly  that  for  most  accurate  solutions  there 
is  an  appropriate  combination  of  interior  stresses  and  boundary  displacements.  For  example, 
in  this  case,  the  linear-moment  and  linear-normal- slope  combination  appears  to  be  a  good 
choice.  Then  the  use  of  more  boundary  displacement  modes  with  the  same  linear  moment 
distribution  can  only  yield  more  inaccurate  solutions  corresponding  to  more  flexible  structures. 
In  fact,  as  we  have  pointed  out  earlier,  the  over- abundance  of  the  boundary  displacement 
modes  will  cause  the  structure  to  be  kinematically  unstable.  It  is  seen,  however,  that  if  the 
approximations  for  the  interior  stresses  and  boundary  displacements  are  improved  simul¬ 
taneously,  the  accuracy  of  the  solution  can  be  improved. 

3.  The  third  example  concerns  a  simply- supported  square  plate  under  uniform  loading. 
For  the  hybrid  model,  a  quadratic  moment  distribution  is  assumed,  hence  unique  generalized 
nodal  forces  can  be  obtained.  The  center  deflection  of  the  plate  has  been  calculated  by  such 
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a  method  using  various  mesh  sizes  as  shown  in  Figure  6.  It  is  seen  that  the  convergence  of 
the  solution  is  extremely  rapid  in  the  uniform  loading  problem.  Solutions  have  also  been 
obtained  by  using  inconsistent  nodal  forces  which  involve  only  vertical  forces  without  moments. 
The  results  are  much  inferior  to  those  obtained  by  using  the  consistently  lumped  nodal  forces. 
Also  included  in  the  figure  for  comparison  are  the  results  by  two  compatible  models  using 
triangular  elements  obtained  by  Clough  and  Tocher  (Reference  9)  and  by  Bazeley,  et  al. 
(Reference  13),  respectively,  and  by  the  equilibrium  model  obtained  by  Morley  (Reference  14). 
In  Morley’ s  formulation,  the  unknown  stress  resultants  (values  of  stress  functions)  are  used 
as  unknowns  and  the  deflection  of  the  plate  can  only  be  calculated  by  integrating  the  moment 
curvature  relations.  In  view  of  the  approximate  character  of  the  moments  provided  by  the 
finite  element  analysis,  the  deflection  obtained  by  integration  is,  in  general,  dependent  upon 
the  chosen  integration  path;  hence  is  not  a  unique  solution.  Thus,  two  different  solutions  by 
Morley  are  shown  in  Figure  6. 

Figure  7  is  a  comparison  of  the  stress  distributions  (M^  along  the  axis  of  symmetry 
y  =  0)  for  this  plate.  It  is  seen  that  the  assumed  stress  methods  (equilibrium  and  hybrid 
models)  can  provide  very  accurate  stress  distributions  while  for  the  compatible  model  only 
stresses  at  the  mid-points  of  the  individual  segments  are  reasonably  close  to  the  exact 
solutions. 
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Figure  6.  Center  Deflection  of  Simply  Supported  Square  Plate  Under  Uniform  Loading 
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SECTION  VI 

TRANSVERSE  SHEAR  EFFECTS  IN  PLATE  BENDING 


The  inclusion  of  transverse  shear  effect  in  the  plate-bending  problems  has  received  only 
limited  attention  in  the  finite  element  application.  Herrmann  (Reference  15)  has  included  such 
effect  in  his  formulation  by  the  Reissner’s  principle.  Severn  and  Taylor  (Reference  11)  have 
pointed  out  the  convenience  of  taking  transverse  shear  effect  into  account  by  the  assumed 
stress  hybrid  model.  Lundgren  (Reference  16)  has  solved  the  problem  of  buckling  of  sandwich 
plates,  the  stiffness  matrices  of  which  are  obtained  by  the  assumed  stress  hybrid  model. 

Only  until  very  recently  was  the  “moderately  thick”  plate-bending  problem  solved  by  the 
finite-element  method  using  compatible  displacement  functions  (Reference  17).  In  that  analysis 
Smith  used  the  plate  theory  advanced  by  Love  so  that  he  can  express  the  strain  energy  in 
terms  of  only  w  and  the  derivatives  of  w.  However,  he  has  to  use  12  degrees  of  freedom  at 
each  node  in  his  rectangular  elements  in  order  to  account  for  the  higher  order  derivatives  of 
w  in  the  strain  energy  functional.  An  alternative  assumed-displacement  approach  for  the 
“moderately  thick”  plate  is  to  distinguish  the  rotations  6X  and  Q  from  the  derivatives  w 
w,y’  resPectively.  use  five  generalized  coordinates  at  each  node.  If,  in  such  a  case, 
only  w,  Qx  and  Q ^  are  used  as  generalized  coordinates  at  the  nodes,  the  displacement  func¬ 
tions  w,  ©x  and  can  only  be  assumed  linear  along  the  boundary  in  order  to  maintain  the 
compatibility  and  if  the  simplest  interpolation  functions  (linear  function  for  a  triangular 
element  and  bilinear  function  for  a  rectangular  element)  are  used,  the  strain  representation 
will  be  extremely  crude.  In  fact,  it  can  be  shown  that  by  such  a  formulation,  the  shear  strain 
energy  will  be  predominating  unless  the  element  size  is  kept  as  small  as  the  plate  thickness. 
Thus  a  finite-element  method  formulated  in  this  manner  will,  in  general,  yield  very  erroneous 
results.  To  improve  the  situation,  the  displacement  functions  must  be  much  more  complicated; 
hence  many  internal  nodes  are  needed. 

When  the  assumed  stress  hybrid  model  is  used,  the  transverse  shear  effect  can  be  taken 
into  account  accurately  by  including  an  adequate  number  of  stress  modes  and  the  use  of  only 
w,  6 x  and  0^  as  generalized  coordinates  can  yield  very  reasonable  solutions.  Figure  8 
indicates  the  results  of  a  finite-element  solution  for  the  center  deflection  of  a  uniformly- 

2 

loaded  sandwich  plate  with  the  parameter  — —  —  4,  where  G  =  shear  modulus  of  the  core, 
h  =  thickness  of  the  core  and  D  =  flexure  modulus  of  the  plate  =  Etfh2/2(1-  v2).  In  this  assumed 
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stress  hybrid  formulation,  the  moment  distributions  are  assumed  quadratic  while  the  boundary 
displacement  w,  8  x  and  9y  are  all  linear.  It  is  seen  that  by  increasing  the  number  of  elements 
the  solution  converges  very  rapidly  to  the  exact  solution  {Reference  18). 


In  the  formulation  of  plate-bending  with  transverse  shear  effect  by  Severn  and  Taylor 
(Reference  11)  and  by  Lundgren  (Reference  16),  the  nodal  displacements  were  still  chosen 
as  w,  w,x  and  w,  with  w  varying  cubically  and  w>n,  linearly  along  each  edge.  The  correct 
choice  should  be^w,  9X,  and  8y  where  9X  and  Qy  should  be  independent  of  w.  Solutions  ob¬ 
tained  by  using  w,  w,x  and  w,y  as  generalized  coordinates  are  also  plotted  in  Figure  8.  It  is 
seen  although  such  a  solution  is  acceptable  for  large  element  sizes  it  becomes  progressively 
inaccurate  when  the  element  size  is  reduced. 


SECTION  VII 
CONCLUSIONS 

The  following  conclusions  can  be  drawn  for  the  assumed  stress  hybrid  method  in  the 
finite-element  analysis: 

(1)  The  method  is  based  on  a  rigorous  variational  principle. 

(2)  This  method  yields  a  structure  more  flexible  than  that  of  a  compatible  model,  but 
more  rigid  than  that  of  an  equilibrium  model,  in  general. 

(3)  A  rational  way  has  been  obtained  for  consistent  lumping  of  body  and  boundary  forces. 

(4)  To  improve  an  element  stiffness  matrix,  the  approximations  of  both  the  interior  stress 
distribution  and  boundary  displacement  distribution  should  be  simultaneously  and  compatibly 
improved. 

(5)  The  method  is  particularly  suitable  to  treat  plate  and  shell  problems  taking  trans¬ 
verse  shear  into  account. 
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Figure  8.  Center^Deflectton  of  Simply-Supported  Square  Sandwich  Plate 


Under  Uniform 
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